Import các thư viện cần thiết
import numpy as np
import pandas as pd
import seaborn as sns
import timeit
import matplotlib.pyplot as plt
import scipy.stats as stats
from sklearn.mixture import GaussianMixture
from sklearn import datasets
from sklearn.cluster import KMeans
from matplotlib.patches import Ellipse
import warnings
warnings.filterwarnings("ignore")
colors = ['purple', 'k', 'r', 'g', 'b']
Ta sẽ implement hàm để tính giá trị này: \begin{equation} \large p(\mathbf x | \mathbf\mu, \mathbf\Sigma) = \frac 1 {({2\pi})^{D/2}|\Sigma|^{1/2}}\exp\left(-\frac 1 2 (\mathbf x -\mathbf\mu)^T\mathbf\Sigma^{-1}(\mathbf x -\mathbf\mu)\right) \end{equation}
# X: (n x d), mu: (d x 1), cov: (d x d) ~ n data points, d dimenions
# return (n x 1), mỗi phần tử ứng với giá trị của hàm trên với n data points
def gaussian(X, mu, cov):
D = X.shape[1]
ret = np.zeros((X.shape[0], 1), dtype = np.float64)
for i in range(X.shape[0]):
tmp = (X[i, :] - mu).T
ret[i] = 1 / ((2 * np.pi) ** (D / 2) * np.linalg.det(cov) ** 0.5) * np.exp(-0.5 * (tmp.T @ np.linalg.inv(cov)) @ tmp)
return ret.reshape(-1, 1)
Ta sẽ initialize các tham số bắt đầu, gồm $\mu$ (dùng K-means), covariance matrix (ma trận đơn vị) và mixing probability $\pi$ (1 / số clusters) của các phân phối Gaussian.
def initialize_param(X, num_clusters):
clusters = []
kmeans = KMeans(num_clusters).fit(X)
mu_k = kmeans.cluster_centers_
for i in range(num_clusters):
clusters.append({
'pi_k': 1.0 / num_clusters,
'mu_k': mu_k[i],
'cov_k': np.identity(X.shape[1], dtype=np.float64)
})
return clusters
Ta sẽ implement thuật toán EM cho GMM, đầu tiên là E-step. Trong bước này ta cần tính $\gamma{(z_{nk})}$ ~ xác suất data point $x_n$ thuộc phân phối Gaussian $k$, khi cho trước $x_n$. \begin{equation} \large \gamma{(z_{nk})}=\frac {\pi_k\mathcal N(\mathbf x_n| \mathbf\mu_k, \mathbf\Sigma_k)}{\sum_{j=1}^K\pi_j\mathcal N(\mathbf x_n| \mathbf\mu_j, \mathbf\Sigma_j)} \end{equation}
def expectation_step(X, clusters):
totals = np.zeros((X.shape[0], 1), dtype=np.float64)
for cluster in clusters:
phi_k = cluster['pi_k']
mu_k = cluster['mu_k']
cov_k = cluster['cov_k']
gamma_nk = (phi_k * gaussian(X, mu_k, cov_k))
totals += gamma_nk
cluster['gamma_nk'] = gamma_nk
for cluster in clusters:
cluster['gamma_nk'] /= totals
return totals
Ta sẽ tới M-step của thuật toán, tìm ra các tham số tối ưu để tối đa hóa hàm log-likelihood. \begin{equation} \large N_k=\sum_{n=1}^N\gamma({z_{nk}}) \end{equation}
\begin{equation} \large \pi_k^{new}=\frac {N_k} N \end{equation}\begin{equation} \large \mu_k^{new}=\frac 1 {N_k} \sum_{n=1}^N\gamma({z_{nk}})\mathbf x_n \end{equation}\begin{equation} \large \Sigma_k^{new}=\frac 1 {N_k} \sum_{n=1}^N\gamma({z_{nk}})(\mathbf x_n-\mathbf\mu_k)(\mathbf x_n-\mathbf\mu_k)^T \end{equation}def maximization_step(X, clusters):
N = float(X.shape[0])
for cluster in clusters:
gamma_nk = cluster['gamma_nk']
cov_k = np.zeros((X.shape[1], X.shape[1]))
N_k = np.sum(gamma_nk, axis=0)
pi_k = N_k / N
mu_k = cluster['mu_k']
for j in range(X.shape[0]):
tmp = (X[j] - mu_k).reshape(-1, 1)
cov_k += gamma_nk[j] * (tmp @ tmp.T)
cov_k /= N_k
mu_k = np.sum(gamma_nk * X, axis=0) / N_k
cluster['pi_k'] = pi_k
cluster['mu_k'] = mu_k
cluster['cov_k'] = cov_k
Ta sẽ implement hàm log-likelihood để đánh giá kết quả. \begin{equation} \large \ln p(\mathbf X)=\sum_{n=1}^N\ln\sum_{k=1}^K\pi_k\mathcal N(\mathbf x_n|\mu_k,\Sigma_k) \end{equation}
def get_log_likelihood(totals):
return np.sum(np.log(totals))
def draw_ellipse(position, covariance, ax=None, **kwargs):
"""Draw an ellipse with a given position and covariance"""
ax = ax or plt.gca()
# Convert covariance to principal axes
if covariance.shape == (2, 2):
U, s, Vt = np.linalg.svd(covariance)
angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
width, height = 2 * np.sqrt(s)
else:
angle = 0
width, height = 2 * np.sqrt(covariance)
# Draw the Ellipse
for nsig in range(1, 4):
ax.add_patch(Ellipse(position, nsig * width, nsig * height,
angle, **kwargs))
def plot_gmm(X, labels, means_, covariances_, weights_, implement=False, label=True, ax=None):
sns.set_style('whitegrid')
ax = ax or plt.gca()
if label:
ax.scatter(X[:, 0], X[:, 1], c=labels, zorder=10, cmap=matplotlib.colors.ListedColormap(colors))
else:
ax.scatter(X[:, 0], X[:, 1], zorder=10)
w_factor = 0.2 / float(max(weights_))
for pos, covar, w in zip(means_, covariances_, weights_):
draw_ellipse(pos, covar, alpha=float(w) * w_factor, color = 'k')
if implement:
plt.axis('square')
plt.show()
def train_gmm(X, num_clusters, num_epochs, threshold, steps=-1):
clusters = initialize_param(X, num_clusters)
likelihoods = np.zeros((num_epochs, ))
predict_prob = np.zeros((X.shape[0], num_clusters))
num_steps = 0
for i in range(num_epochs):
totals = expectation_step(X, clusters)
if num_steps == steps or (i == 0 and steps != -1):
num_steps = 0
means_ = []
covariances_ = []
weights_ = []
for j, cluster in enumerate(clusters):
predict_prob[:, j] = cluster['gamma_nk'].reshape(-1)
means_.append(cluster['mu_k'])
covariances_.append(cluster['cov_k'])
weights_.append(cluster['pi_k'])
labels = np.zeros((X.shape[0], ))
for j in range(X.shape[0]):
labels[j] = np.argmax(predict_prob[j])
plot_gmm(X, labels, means_, covariances_, weights_, 1)
maximization_step(X, clusters)
num_steps += 1
likelihoods[i] = get_log_likelihood(totals)
if abs(likelihoods[i] - likelihoods[i - 1]) < threshold:
print('Improvement smaller than threshold! Finished Training')
likelihoods = np.resize(likelihoods, i + 1)
break
print('Epoch: ', i + 1, 'Likelihood: ', likelihoods[i])
for i, cluster in enumerate(clusters):
predict_prob[:, i] = cluster['gamma_nk'].reshape(-1)
return clusters, likelihoods, predict_prob
def plot_log_likelihood(likelihoods):
sns.set_style('whitegrid')
plt.figure(figsize=(5, 5))
plt.title('Log-Likelihood')
plt.plot(np.arange(1, len(likelihoods) + 1), likelihoods)
plt.show()
Ta sẽ tạo ra một dataset 1D bằng một số phân phối Gaussian, và dùng GMM để xác định lại các phân phối Gaussian đã tạo ra model đó.
# Generate 1D data using many Gaussian distributions
def gen_1D_data_with_gds(mu, std, sz):
num_gd = len(mu)
data = np.random.normal(mu[0], std[0], sz[0])
for i in range(1, num_gd):
data = np.concatenate((data, np.random.normal(mu[i], std[i], sz[i])), axis = 0)
return data
mu = [0, 1, 3]
std = [0.35, 0.2, 0.3]
sz = [2000, 2000, 3000]
data_1D_2gd = gen_1D_data_with_gds(mu, std, sz).reshape(-1, 1)
sns.set_style('whitegrid')
plt.hist(data_1D_2gd, bins=100, density=1, color="lightblue")
plt.title('1D dataset generated by 3 Gaussian distributions')
plt.show()
Ta sẽ so sánh khi sử dụng GMM của sklearn và GMM self-implementation cho dataset 1D đã tạo ra ở trên.
print('Sử dụng GMM của thư viện sklearn')
gmm = GaussianMixture(n_components=3, verbose=2, random_state=42)
start = timeit.default_timer()
gmm.fit(data_1D_2gd)
stop = timeit.default_timer()
print('Time of GMM sklearn:', stop - start)
sns.set_style('whitegrid')
plt.hist(data_1D_2gd, bins=100, density=1, color="lightblue", zorder=0)
x_axis = np.linspace(-1.5, 4.0, 100)
for i in range(3):
plt.plot(x_axis, gmm.weights_[i] * stats.norm.pdf(x_axis, gmm.means_[i, 0], np.sqrt(gmm.covariances_[i, 0])[0]), linestyle='solid')
print("Gaussian distribution {}: Mean: {}, Std: {}".format(i + 1, gmm.means_[i, 0], np.sqrt(gmm.covariances_[i, 0])[0]))
plt.show()
Sử dụng GMM của thư viện sklearn Initialization 0 Initialization converged: True time lapse 0.02600s ll -1.17533 Time of GMM sklearn: 0.0305377 Gaussian distribution 1: Mean: -0.039891987741650994, Std: 0.3308096901049394 Gaussian distribution 2: Mean: 3.0083572742343283, Std: 0.29778797934389883 Gaussian distribution 3: Mean: 0.9758018123747146, Std: 0.21922464533712888
print('Sử dụng GMM self-implementation')
num_clusters = 3
num_epochs = 100
threshold = 1e-9
start = timeit.default_timer()
clusters, likelihoods, predict_prob = train_gmm(data_1D_2gd, num_clusters, num_epochs, threshold)
stop = timeit.default_timer()
print('Time of GMM self-implementation:', stop - start)
plt.hist(data_1D_2gd, bins=100, density=1, color="lightblue")
x_axis = np.linspace(-1.5, 4.0, 100)
for i in range(num_clusters):
plt.plot(x_axis, clusters[i]['pi_k'] * stats.norm.pdf(x_axis, clusters[i]['mu_k'][0], np.sqrt(clusters[i]['cov_k'][0, 0])), linestyle='solid')
print("Gaussian distribution {}: Mean: {}, Std: {}".format(i + 1, clusters[i]['mu_k'][0], np.sqrt(clusters[i]['cov_k'][0, 0])))
plt.show()
plot_log_likelihood(likelihoods)
Sử dụng GMM self-implementation Epoch: 1 Likelihood: -11879.900145255608 Epoch: 2 Likelihood: -10415.54440492494 Epoch: 3 Likelihood: -9209.42034115552 Epoch: 4 Likelihood: -8941.08072592746 Epoch: 5 Likelihood: -8886.67623982393 Epoch: 6 Likelihood: -8869.113939127197 Epoch: 7 Likelihood: -8864.223093157292 Epoch: 8 Likelihood: -8862.547048760087 Epoch: 9 Likelihood: -8861.532025099457 Epoch: 10 Likelihood: -8860.532196663487 Epoch: 11 Likelihood: -8859.351209753244 Epoch: 12 Likelihood: -8857.88246226405 Epoch: 13 Likelihood: -8856.014262790366 Epoch: 14 Likelihood: -8853.59266190622 Epoch: 15 Likelihood: -8850.387811628505 Epoch: 16 Likelihood: -8846.042624842972 Epoch: 17 Likelihood: -8839.982583569246 Epoch: 18 Likelihood: -8831.250080749305 Epoch: 19 Likelihood: -8818.193555906553 Epoch: 20 Likelihood: -8797.88686850547 Epoch: 21 Likelihood: -8765.119218438176 Epoch: 22 Likelihood: -8711.09917210203 Epoch: 23 Likelihood: -8624.027614903735 Epoch: 24 Likelihood: -8499.702503326265 Epoch: 25 Likelihood: -8367.623941490301 Epoch: 26 Likelihood: -8281.7241019684 Epoch: 27 Likelihood: -8247.649329322896 Epoch: 28 Likelihood: -8235.840240802503 Epoch: 29 Likelihood: -8230.207683172483 Epoch: 30 Likelihood: -8226.52974446038 Epoch: 31 Likelihood: -8223.91024351924 Epoch: 32 Likelihood: -8222.053886360036 Epoch: 33 Likelihood: -8220.766675118864 Epoch: 34 Likelihood: -8219.892976264706 Epoch: 35 Likelihood: -8219.310532814308 Epoch: 36 Likelihood: -8218.927858714153 Epoch: 37 Likelihood: -8218.679322390466 Epoch: 38 Likelihood: -8218.51936553813 Epoch: 39 Likelihood: -8218.417149266763 Epoch: 40 Likelihood: -8218.352194148722 Epoch: 41 Likelihood: -8218.311097212727 Epoch: 42 Likelihood: -8218.285184091252 Epoch: 43 Likelihood: -8218.26888867127 Epoch: 44 Likelihood: -8218.25866285646 Epoch: 45 Likelihood: -8218.252256470912 Epoch: 46 Likelihood: -8218.248248130316 Epoch: 47 Likelihood: -8218.24574275414 Epoch: 48 Likelihood: -8218.24417804903 Epoch: 49 Likelihood: -8218.243201447352 Epoch: 50 Likelihood: -8218.242592210534 Epoch: 51 Likelihood: -8218.242212297324 Epoch: 52 Likelihood: -8218.241975460995 Epoch: 53 Likelihood: -8218.241827854199 Epoch: 54 Likelihood: -8218.241735876854 Epoch: 55 Likelihood: -8218.241678572243 Epoch: 56 Likelihood: -8218.241642874045 Epoch: 57 Likelihood: -8218.241620637771 Epoch: 58 Likelihood: -8218.241606787908 Epoch: 59 Likelihood: -8218.241598162027 Epoch: 60 Likelihood: -8218.241592789964 Epoch: 61 Likelihood: -8218.241589444446 Epoch: 62 Likelihood: -8218.241587361046 Epoch: 63 Likelihood: -8218.24158606365 Epoch: 64 Likelihood: -8218.241585255737 Epoch: 65 Likelihood: -8218.241584752644 Epoch: 66 Likelihood: -8218.241584439367 Epoch: 67 Likelihood: -8218.241584244288 Epoch: 68 Likelihood: -8218.241584122816 Epoch: 69 Likelihood: -8218.241584047173 Epoch: 70 Likelihood: -8218.241584000076 Epoch: 71 Likelihood: -8218.241583970746 Epoch: 72 Likelihood: -8218.241583952484 Epoch: 73 Likelihood: -8218.241583941111 Epoch: 74 Likelihood: -8218.241583934032 Epoch: 75 Likelihood: -8218.241583929623 Epoch: 76 Likelihood: -8218.241583926878 Epoch: 77 Likelihood: -8218.241583925166 Epoch: 78 Likelihood: -8218.241583924102 Improvement smaller than threshold! Finished Training Time of GMM self-implementation: 37.9598316 Gaussian distribution 1: Mean: 3.0083405660957885, Std: 0.2978173251358502 Gaussian distribution 2: Mean: -0.0141648827590889, Std: 0.3501462810105499 Gaussian distribution 3: Mean: 0.9894713552543288, Std: 0.20839597070036683
Ta thấy cả 2 GMM đều tìm được 2 phân phối Gaussian khá chính xác so với 2 phân phối ban đầu dùng để sinh ra data. Mặc dù vậy GMM của sklearn chạy với tốc độ nhanh gấp nhiều lần trong khi GMM khi tự implement lại thì chạy khá chậm.
Giờ ta hãy nghịch thử với phân phối điểm thi thptqg năm 2021 (xin gửi lời cảm ơn đến anh mentor VLTAnh vì data ạ :>)
# dt21 = pd.read_csv('thptqg2021', sep = '\t')
# %store dt21
%store -r dt21
dt21.head()
dt21.columns.values[4:]
array(['Toán', 'Văn', 'Lý', 'Hoá', 'Sinh', 'KHTN', 'Lịch Sử', 'Địa Lý',
'GDCD', 'KHXH', 'Ngoại Ngữ'], dtype=object)
subjects = dt21.columns.values[4:]
i = 0;
for sub in subjects:
i += 1
sns.set_style('whitegrid')
sns.distplot(pd.to_numeric(dt21[sub], errors='coerce'), bins = 40, color='lightblue')
plt.show()
Ta sẽ thử áp dụng GMM cho phân phối điểm môn Ngoại Ngữ (vì là môn duy nhất có 2 đỉnh)
print('Sử dụng GMM của thư viện sklearn')
dtNN = pd.to_numeric(dt21['Ngoại Ngữ'], errors='coerce').values;
dtNN = dtNN[~np.isnan(dtNN)].reshape(-1, 1)
gmm = GaussianMixture(n_components=2, verbose=2, random_state=42)
start = timeit.default_timer()
gmm.fit(dtNN)
stop = timeit.default_timer()
print('Time of GMM sklearn:', stop - start)
sns.set_style('whitegrid')
plt.hist(dtNN, bins=40, density=1, color="lightblue")
x_axis = np.linspace(-0.3, 11, 40)
for i in range(2):
plt.plot(x_axis, gmm.weights_[i] * stats.norm.pdf(x_axis, gmm.means_[i, 0], np.sqrt(gmm.covariances_[i, 0])))
print("Gaussian distribution {}: Mean: {}, Std: {}".format(i + 1, gmm.means_[i, 0], np.sqrt(gmm.covariances_[i, 0])[0]))
plt.title('Ngoại Ngữ', y=-0.16)
plt.show()
Sử dụng GMM của thư viện sklearn Initialization 0 Initialization converged: True time lapse 0.79669s ll -2.12194 Time of GMM sklearn: 0.9174209999999903 Gaussian distribution 1: Mean: 4.263676305632633, Std: 1.197547321924117 Gaussian distribution 2: Mean: 8.034974117926646, Std: 1.1470684546266454
print('Sử dụng GMM self-implementation')
num_clusters = 2
num_epochs = 30
threshold = 1e-9
start = timeit.default_timer()
clusters, likelihoods, predict_prob = train_gmm(dtNN, num_clusters, num_epochs, threshold)
stop = timeit.default_timer()
print('Time of GMM self-implementation:', stop - start)
sns.set_style('whitegrid')
plt.hist(dtNN, bins=40, density=1, color="lightblue")
x_axis = np.linspace(-0.3, 11.0, 40)
for i in range(num_clusters):
plt.plot(x_axis, clusters[i]['pi_k'] * stats.norm.pdf(x_axis, clusters[i]['mu_k'][0], np.sqrt(clusters[i]['cov_k'][0, 0])), linestyle='solid')
print("Gaussian distribution {}: Mean: {}, Std: {}".format(i + 1, clusters[i]['mu_k'][0], np.sqrt(clusters[i]['cov_k'][0, 0])))
plt.show()
plot_log_likelihood(likelihoods)
Sử dụng GMM self-implementation Epoch: 1 Likelihood: -1746951.2355478865 Epoch: 2 Likelihood: -1725959.349239038 Epoch: 3 Likelihood: -1724710.6990825243 Epoch: 4 Likelihood: -1724382.2927940225 Epoch: 5 Likelihood: -1724159.43939155 Epoch: 6 Likelihood: -1723957.2409362076 Epoch: 7 Likelihood: -1723765.1104538152 Epoch: 8 Likelihood: -1723581.555423499 Epoch: 9 Likelihood: -1723406.2078877126 Epoch: 10 Likelihood: -1723238.817073564 Epoch: 11 Likelihood: -1723079.134125698 Epoch: 12 Likelihood: -1722926.9050382355 Epoch: 13 Likelihood: -1722781.8734092826 Epoch: 14 Likelihood: -1722643.7826795264 Epoch: 15 Likelihood: -1722512.3775804087 Epoch: 16 Likelihood: -1722387.4051645258 Epoch: 17 Likelihood: -1722268.615617903 Epoch: 18 Likelihood: -1722155.7629370904 Epoch: 19 Likelihood: -1722048.6055053142 Epoch: 20 Likelihood: -1721946.906583757 Epoch: 21 Likelihood: -1721850.4347276092 Epoch: 22 Likelihood: -1721758.964133557 Epoch: 23 Likelihood: -1721672.274925081 Epoch: 24 Likelihood: -1721590.1533807588 Epoch: 25 Likelihood: -1721512.3921108702 Epoch: 26 Likelihood: -1721438.7901871132 Epoch: 27 Likelihood: -1721369.1532299917 Epoch: 28 Likelihood: -1721303.2934581602 Epoch: 29 Likelihood: -1721241.0297037815 Epoch: 30 Likelihood: -1721182.187397595 Time of GMM self-implementation: 1109.1125655 Gaussian distribution 1: Mean: 8.288647265813061, Std: 0.9837302392335664 Gaussian distribution 2: Mean: 4.473478953087329, Std: 1.3276176572024334
Vì số lượng lên tới hơn 800000 học sinh, GMM self-implementation chạy chậm hơn nhiều so với thư viện sklearn.
Ta sẽ viết 2 hàm dùng cho datasets 2D lần lượt là GMM của sklearn và GMM self-implementation
def gmm_sklearn_2D_dataset(X, num_clusters, label=False, square=1):
gmm = GaussianMixture(n_components=num_clusters, verbose=2, random_state=42)
start = timeit.default_timer()
gmm.fit(X)
stop = timeit.default_timer()
print('Time of GMM sklearn:', stop - start)
labels = gmm.predict(X)
plot_gmm(X, labels, gmm.means_, gmm.covariances_, gmm.weights_, square, label)
def gmm_2D_dataset(X, num_clusters, steps=30, num_epochs=100, threshold=1e-9):
start = timeit.default_timer()
clusters, likelihoods, predict_prob = train_gmm(X, num_clusters, num_epochs, threshold, steps)
stop = timeit.default_timer()
print('Time of GMM self-implementation:', stop - start)
labels = np.zeros((X.shape[0], ))
for i in range(X.shape[0]):
labels[i] = np.argmax(predict_prob[i])
means_ = []
covariances_ = []
weights_ = []
for cluster in clusters:
means_.append(cluster['mu_k'])
covariances_.append(cluster['cov_k'])
weights_.append(cluster['pi_k'])
plot_gmm(X, labels, means_, covariances_, weights_, 1)
plot_log_likelihood(likelihoods)
Ta sẽ dùng sklearn.datasets để tạo ra một data blobs với 1000 điểm và transform các điểm thành data như bên dưới.
X, y = datasets.make_blobs(n_samples=1000, random_state=170)
transformation = [[0.6, -0.5], [-0.5, 0.9]]
X_aniso = np.dot(X, transformation)
sns.set_style('whitegrid')
plt.scatter(X_aniso[:, 0], X_aniso[:, 1], c = y, cmap=matplotlib.colors.ListedColormap(colors))
plt.axis('square')
plt.show()
gmm_sklearn_2D_dataset(X_aniso, 3, 1)
Initialization 0 Initialization converged: True time lapse 0.00000s ll -2.63988 Time of GMM sklearn: 0.011501799999905415
gmm_2D_dataset(X_aniso, 3, 1)
# sau 3 epochs thì plot các Gaussian đã tìm được ở bước hiện tại
Epoch: 1 Likelihood: -3498.697459025682
Epoch: 2 Likelihood: -3163.959055888265
Epoch: 3 Likelihood: -2992.7946958109596
Epoch: 4 Likelihood: -2884.525903829919
Epoch: 5 Likelihood: -2692.727067598239
Epoch: 6 Likelihood: -2640.3812948562704
Epoch: 7 Likelihood: -2639.8859323768293
Epoch: 8 Likelihood: -2639.8810061437935
Epoch: 9 Likelihood: -2639.8809021072707
Epoch: 10 Likelihood: -2639.8808998876075
Epoch: 11 Likelihood: -2639.8808998402355
Epoch: 12 Likelihood: -2639.8808998392233
Improvement smaller than threshold! Finished Training Time of GMM self-implementation: 3.109119899999996
Ta thấy cả 2 GMM đều đã phân loại được đúng datasets, điều mà thuật toán clustering như K-means phân loại không chính xác.
Nhìn các Gaussian sau các bước của thuật toán EM thay đổi, ta thấy rõ cách mà GMM cải tiến hơn so với K-means.
Tiếp theo, ta sẽ chạy trên một dataset 2D mà các thuật toán như K-means cũng không phân loại một cách chính xác được.
X_blobs, y = datasets.make_blobs(n_samples=3000,
cluster_std=[1.0, 2.5, 0.5],
random_state=170)
sns.set_style('whitegrid')
plt.scatter(X_blobs[:, 0], X_blobs[:, 1], c = y, cmap=matplotlib.colors.ListedColormap(colors))
plt.axis('square')
plt.show()
gmm_sklearn_2D_dataset(X_blobs, 3, 1)
Initialization 0 Initialization converged: True time lapse 0.01562s ll -4.00599 Time of GMM sklearn: 0.02225670000007085
gmm_2D_dataset(X_blobs, 3, 5)
Epoch: 1 Likelihood: -14490.116869862377 Epoch: 2 Likelihood: -12606.617088183435 Epoch: 3 Likelihood: -12134.789144176633 Epoch: 4 Likelihood: -12033.261945451317 Epoch: 5 Likelihood: -12019.96823729699
Epoch: 6 Likelihood: -12018.011941221368 Epoch: 7 Likelihood: -12017.599867156205 Epoch: 8 Likelihood: -12017.492624044311 Epoch: 9 Likelihood: -12017.4623292978 Epoch: 10 Likelihood: -12017.453509762167
Epoch: 11 Likelihood: -12017.450911984546 Epoch: 12 Likelihood: -12017.450143014585 Epoch: 13 Likelihood: -12017.449914870764 Epoch: 14 Likelihood: -12017.449847107193 Epoch: 15 Likelihood: -12017.449826968339
Epoch: 16 Likelihood: -12017.449820981385 Epoch: 17 Likelihood: -12017.449819201272 Epoch: 18 Likelihood: -12017.449818671941 Epoch: 19 Likelihood: -12017.449818514533 Epoch: 20 Likelihood: -12017.449818467723
Epoch: 21 Likelihood: -12017.449818453802 Epoch: 22 Likelihood: -12017.44981844966 Epoch: 23 Likelihood: -12017.449818448431 Improvement smaller than threshold! Finished Training Time of GMM self-implementation: 6.034029700000019
GMM nhóm cụm một cách chính xác mặc dù phân phối của 3 cụm có độ lệch chuẩn khác nhau rõ rệt.
Tiếp đến ta sẽ xem một số datasets 2D mà GMM chưa thể phân cụm chính xác
X_moon, y_moon = datasets.make_moons(500, noise=0.05, random_state=42)
sns.set_style('whitegrid')
plt.scatter(X_moon[:, 0], X_moon[:, 1], c=y_moon, cmap=matplotlib.colors.ListedColormap(colors));
gmm_sklearn_2D_dataset(X_moon, 2, 1, 0)
Initialization 0 Initialization converged: True time lapse 0.00699s ll -1.71401 Time of GMM sklearn: 0.007359899999983099
X_cir, y_cir = datasets.make_circles(n_samples=1500, factor=.5, noise=.05)
sns.set_style('whitegrid')
plt.scatter(X_cir[:, 0], X_cir[:, 1], c=y_cir, cmap=matplotlib.colors.ListedColormap(colors));
plt.axis('square')
plt.show()
gmm_sklearn_2D_dataset(X_cir, 2, 1)
Initialization 0 Initialization converged: True time lapse 0.01138s ll -1.64053 Time of GMM sklearn: 0.00998200000003635
Ta dễ dàng thấy đối với các datasets mà các điểm dữ liệu có phân phối đặc biệt (không tuân theo phân phối Gaussian) thì GMM không thể phân cụm chính xác.
Ngoài ứng dụng trong các bài toán clustering, GMM còn dược sử dụng cho density estimation nhằm mô tả phân phối của data, từ đó được dùng như một generative model giúp sinh các data mới theo phân phối của data cho trước.
Ta hãy cùng quay lại với datasets moon ở trên.
X_moon, y_moon = datasets.make_moons(500, noise=0.05, random_state=42)
sns.set_style('whitegrid')
plt.scatter(X_moon[:, 0], X_moon[:, 1], c=y_moon, cmap=matplotlib.colors.ListedColormap(colors));
Khi ta dùng GMM với số lượng cụm khoảng 8, thì ta sẽ được một mô tả khá chính xác phân phối của dữ liệu, ta sẽ dùng các phân phối Gaussian đó để sinh ra các dữ liệu mới.
gmm = GaussianMixture(n_components=8, verbose=2, random_state=42)
gmm.fit(X_moon)
sns.set_style('whitegrid')
plot_gmm(X_moon, [], gmm.means_, gmm.covariances_, gmm.weights_, 0, 0)
Initialization 0 Iteration 10 time lapse 0.01562s ll change 0.00155 Initialization converged: True time lapse 0.01562s ll -0.42189
# Generate các data mới theo phân phối của data trên
X_new, y_new = gmm.sample(500)
sns.set_style('whitegrid')
plt.scatter(X_new[:, 0], X_new[:, 1]);
500 điểm data mới được sinh ra dựa trên 8 Gaussian model mà GMM đã phân cụm.
Ta sẽ dễ dàng thấy rằng nếu ta tăng số cụm của GMM lên càng nhiều, thì log-likelihood sẽ càng tăng. Nhưng đôi khi việc phân bằng số cụm nhiều thì model của ta sẽ không đủ "general" để đưa ra được các cách nhóm cụm hữu ích, hoặc nếu số cụm ít quá thì model sẽ không đủ "detailed" để mô tả các đặc điểm nổi bật của data. Vì vậy, ta sẽ sử dụng AIC và BIC để có thể chọn số cụm sao cho phù hợp với data.
models = [GaussianMixture(num, random_state=42).fit(X_moon)
for num in range(1, 21)]
sns.set_style('whitegrid')
plt.plot([m.bic(X_moon) for m in models], label='BIC')
plt.plot([m.aic(X_moon) for m in models], label='AIC')
plt.xticks(range(0,21, 4))
plt.legend()
plt.xlabel('num_clusters');
Ta thấy số cụm nên sử dụng cho GMM để mô tả hiệu quả phân phối X_moon là vào khoảng từ 8 đến 12, vùng mà minimize giá trị AIC và BIC.
Ta sẽ lấy data handwritten digits từ thư viện sklearn, mỗi digit sẽ được biểu diễn bằng 64 pixels. Sau đó ta sẽ sử dụng GMM như một generative model để sinh ra các digits mới.
digits = datasets.load_digits()
digits.data.shape
(1797, 64)
def plot_digits(data):
sns.set_style('whitegrid')
fig, ax = plt.subplots(10, 10, figsize=(8, 8),
subplot_kw=dict(xticks=[], yticks=[]))
fig.subplots_adjust(hspace=0.05, wspace=0.05)
for i, axi in enumerate(ax.flat):
im = axi.imshow(data[i].reshape(8, 8), cmap='binary')
im.set_clim(0, 16)
plot_digits(digits.data)
models = [GaussianMixture(n, random_state=42)
for n in range(30, 210, 10)]
sns.set_style('whitegrid')
plt.plot(range(30, 210, 10), [model.fit(digits.data).aic(digits.data) for model in models], label='AIC');
plt.plot(range(30, 210, 10), [model.fit(digits.data).bic(digits.data) for model in models], label='BIC');
plt.legend()
plt.show()
Ta sẽ sử dụng GMM với 100 cụm
gmm = GaussianMixture(n_components=100, verbose=2, random_state=42)
gmm.fit(digits.data)
Initialization 0 Initialization converged: True time lapse 0.44156s ll 215.66135
GaussianMixture(n_components=100, random_state=42, verbose=2)
new_data, _ = gmm.sample(100)
plot_digits(new_data)
Ta hãy cùng thử nghiệm mô hình thuật toán EM áp dụng cho các bài toán GMM semi-supervised, khi mà có cả labelled và unlabelled examples.
# X: (n x d), mu: (d x 1), cov: (d x d) ~ n data points, d dimenions
# return (n x 1), mỗi phần tử ứng với giá trị của hàm trên với n data points
def gaussian(X, mu, cov):
D = X.shape[1]
ret = np.zeros((X.shape[0], 1), dtype = np.float64)
for i in range(X.shape[0]):
tmp = (X[i, :] - mu).T
ret[i] = 1 / ((2 * np.pi) ** (D / 2) * np.linalg.det(cov) ** 0.5) * np.exp(-0.5 * (tmp.T @ np.linalg.inv(cov)) @ tmp)
return ret.reshape(-1, 1)
def initialize_param(X, X_label, y_label, num_clusters):
clusters = []
kmeans = KMeans(num_clusters).fit(X)
mu_k = kmeans.cluster_centers_
for i in range(num_clusters):
clusters.append({
'pi_k': 1.0 / num_clusters,
'mu_k': mu_k[i],
'cov_k': np.identity(X.shape[1], dtype=np.float64),
'id' : i
})
idx = np.arange(0, y_label.shape[0])
for i in range(np.max(y_label)):
clusters[i]['mu_k'] = X_label[np.random.choice(idx[y_label == i], size = 1)]
return clusters
def expectation_step(X, clusters):
totals = np.zeros((X.shape[0], 1), dtype=np.float64)
for cluster in clusters:
phi_k = cluster['pi_k']
mu_k = cluster['mu_k']
cov_k = cluster['cov_k']
gamma_nk = (phi_k * gaussian(X, mu_k, cov_k))
totals += gamma_nk
cluster['gamma_nk'] = gamma_nk
for cluster in clusters:
cluster['gamma_nk'] /= totals
def maximization_step(X, X_label, y_label, clusters, alpha):
N = float(X.shape[0])
N_label = float(X_label.shape[0])
for cluster in clusters:
gamma_nk = cluster['gamma_nk']
cov_k = np.zeros((X.shape[1], X.shape[1]))
N_k = np.sum(gamma_nk, axis=0) + alpha * np.sum(y_label == cluster['id'])
pi_k = N_k / (N + alpha * N_label)
mu_k = cluster['mu_k']
for j in range(X.shape[0]):
tmp = (X[j] - mu_k).reshape(-1, 1)
cov_k += gamma_nk[j] * (tmp @ tmp.T)
for j in range(X_label.shape[0]):
if y_label[j] == cluster['id']:
tmp = (X_label[j] - mu_k).reshape(-1, 1)
cov_k += alpha * (tmp @ tmp.T)
cov_k /= N_k
mu_k = (np.sum(gamma_nk * X, axis=0) + alpha * np.sum(X_label[y_label == cluster['id']], axis=0)) / N_k
cluster['pi_k'] = pi_k
cluster['mu_k'] = mu_k
cluster['cov_k'] = cov_k
def train_semi_supervised_gmm(X, X_label, y_label, num_clusters, num_epochs, alpha):
clusters = initialize_param(X, X_label, y_label, num_clusters)
predict_prob = np.zeros((X.shape[0], num_clusters))
for i in range(num_epochs):
expectation_step(X, clusters)
maximization_step(X, X_label, y_label, clusters, alpha)
print('Epoch: ', i + 1)
for i, cluster in enumerate(clusters):
predict_prob[:, i] = cluster['gamma_nk'].reshape(-1)
return clusters, predict_prob
Ta sẽ tạo ra một dataset 2D để so sánh giữa Unsupervised GMM và Semi-Supervised GMM.
# Generate 2D data using many Multivariate Gaussian distributions
def gen_2D_data_with_gds(mu, cov, sz):
num_gd = len(mu)
X = np.zeros((np.sum(sz), 2))
y = np.zeros((np.sum(sz), ), dtype = np.int)
lst = 0
for i in range(num_gd):
X[lst : lst + sz[i], :] = np.random.multivariate_normal(mu[i], cov[i], sz[i])
y[lst : lst + sz[i]] = i
lst += sz[i]
return X, y
mu = [[5, 5], [7, 7.5], [8, 6], [7, 5.2], [3, 3]]
cov = [[[2, 0], [1, 2]], [[0.2, 0.1], [0.1, 0.2]], [[0.2, -0.05], [-0.1, 0.2]], [[0.2, -0.2], [0.1, 0.2]], [[0.2, -0.2], [0.1, 0.2]]]
sz = [400, 300, 180, 150, 70]
X, _ = gen_2D_data_with_gds(mu, cov, sz)
sz = [10, 10, 10, 10, 10]
X_label, y_label = gen_2D_data_with_gds(mu, cov, sz)
sns.set_style('whitegrid')
plt.scatter(X[:, 0], X[:, 1])
plt.scatter(X_label[:, 0], X_label[:, 1], c = y_label, cmap = matplotlib.colors.ListedColormap(colors))
plt.axis('square')
plt.show()
Ta thấy dataset 2D ở trên được generated với 5 phân phối Gaussian, mỗi phân phối ta biết tầm 10 điểm thuộc phân phối đó. (Những điểm màu khác màu đỏ)
gmm_sklearn_2D_dataset(X, 5, 1)
Initialization 0 Initialization converged: True time lapse 0.01300s ll -3.22524 Time of GMM sklearn: 0.015185299999984636
Ta thấy GMM với thuật toán EM truyền thống không phân cụm đúng với các phân phối Gaussian mà ta đã dùng để generate dataset này.
clusters, predict_prob = train_semi_supervised_gmm(X, X_label, y_label, 5, 100, 1.2)
labels = np.zeros((X.shape[0], ))
for i in range(X.shape[0]):
labels[i] = np.argmax(predict_prob[i])
means_ = []
covariances_ = []
weights_ = []
for cluster in clusters:
means_.append(cluster['mu_k'])
covariances_.append(cluster['cov_k'])
weights_.append(cluster['pi_k'])
plot_gmm(X, labels, means_, covariances_, weights_, 1)
Epoch: 1 Epoch: 2 Epoch: 3 Epoch: 4 Epoch: 5 Epoch: 6 Epoch: 7 Epoch: 8 Epoch: 9 Epoch: 10 Epoch: 11 Epoch: 12 Epoch: 13 Epoch: 14 Epoch: 15 Epoch: 16 Epoch: 17 Epoch: 18 Epoch: 19 Epoch: 20 Epoch: 21 Epoch: 22 Epoch: 23 Epoch: 24 Epoch: 25 Epoch: 26 Epoch: 27 Epoch: 28 Epoch: 29 Epoch: 30 Epoch: 31 Epoch: 32 Epoch: 33 Epoch: 34 Epoch: 35 Epoch: 36 Epoch: 37 Epoch: 38 Epoch: 39 Epoch: 40 Epoch: 41 Epoch: 42 Epoch: 43 Epoch: 44 Epoch: 45 Epoch: 46 Epoch: 47 Epoch: 48 Epoch: 49 Epoch: 50 Epoch: 51 Epoch: 52 Epoch: 53 Epoch: 54 Epoch: 55 Epoch: 56 Epoch: 57 Epoch: 58 Epoch: 59 Epoch: 60 Epoch: 61 Epoch: 62 Epoch: 63 Epoch: 64 Epoch: 65 Epoch: 66 Epoch: 67 Epoch: 68 Epoch: 69 Epoch: 70 Epoch: 71 Epoch: 72 Epoch: 73 Epoch: 74 Epoch: 75 Epoch: 76 Epoch: 77 Epoch: 78 Epoch: 79 Epoch: 80 Epoch: 81 Epoch: 82 Epoch: 83 Epoch: 84 Epoch: 85 Epoch: 86 Epoch: 87 Epoch: 88 Epoch: 89 Epoch: 90 Epoch: 91 Epoch: 92 Epoch: 93 Epoch: 94 Epoch: 95 Epoch: 96 Epoch: 97 Epoch: 98 Epoch: 99 Epoch: 100
Trong khi, Semi-Supervised GMM phân cụm rất chính xác, chỉ nhờ vào 50 điểm cho trước thuộc các phân phối Gaussian
Ta hãy cùng đến với một dataset khác.
# data = pd.read_csv('semi.csv')
# %store data
%store -r data
data.head()
| x_1 | x_2 | z | |
|---|---|---|---|
| 0 | -0.664166 | 1.051539 | -1 |
| 1 | 0.321250 | 2.966127 | 3 |
| 2 | -0.078748 | 0.955693 | -1 |
| 3 | -0.036538 | 1.030159 | -1 |
| 4 | -1.553956 | 2.776952 | -1 |
X1 = data[data['z'] == -1].drop(['z'], axis = 1).to_numpy()
X1_label = data[data['z'] != -1].drop(['z'], axis = 1).to_numpy()
y1_label = data['z'][data['z'] != -1].to_numpy()
sns.set_style('whitegrid')
plt.scatter(X1[:, 0], X1[:, 1])
plt.scatter(X1_label[:, 0], X1_label[:, 1], c = y1_label, cmap = matplotlib.colors.ListedColormap(colors))
plt.axis('square')
plt.show()
Ta thấy dataset 2D ở trên được generated với 4 phân phối Gaussian, mỗi phân phối ta chỉ biết 5 điểm thuộc phân phối đó. (Những điểm màu khác màu đỏ)
gmm_sklearn_2D_dataset(X1, 4, 1)
Initialization 0 Iteration 10 time lapse 0.01562s ll change 0.00209 Initialization converged: True time lapse 0.01562s ll -1.83616 Time of GMM sklearn: 0.021091500000011365
Ta thấy GMM với thuật toán EM truyền thống không phân cụm đúng với các phân phối Gaussian mà ta đã dùng để generate dataset này.
clusters, predict_prob = train_semi_supervised_gmm(X1, X1_label, y1_label, 4, 100, 1.2)
labels = np.zeros((X1.shape[0], ))
for i in range(X1.shape[0]):
labels[i] = np.argmax(predict_prob[i])
means_ = []
covariances_ = []
weights_ = []
for cluster in clusters:
means_.append(cluster['mu_k'])
covariances_.append(cluster['cov_k'])
weights_.append(cluster['pi_k'])
plot_gmm(X1, labels, means_, covariances_, weights_, 1)
Epoch: 1 Epoch: 2 Epoch: 3 Epoch: 4 Epoch: 5 Epoch: 6 Epoch: 7 Epoch: 8 Epoch: 9 Epoch: 10 Epoch: 11 Epoch: 12 Epoch: 13 Epoch: 14 Epoch: 15 Epoch: 16 Epoch: 17 Epoch: 18 Epoch: 19 Epoch: 20 Epoch: 21 Epoch: 22 Epoch: 23 Epoch: 24 Epoch: 25 Epoch: 26 Epoch: 27 Epoch: 28 Epoch: 29 Epoch: 30 Epoch: 31 Epoch: 32 Epoch: 33 Epoch: 34 Epoch: 35 Epoch: 36 Epoch: 37 Epoch: 38 Epoch: 39 Epoch: 40 Epoch: 41 Epoch: 42 Epoch: 43 Epoch: 44 Epoch: 45 Epoch: 46 Epoch: 47 Epoch: 48 Epoch: 49 Epoch: 50 Epoch: 51 Epoch: 52 Epoch: 53 Epoch: 54 Epoch: 55 Epoch: 56 Epoch: 57 Epoch: 58 Epoch: 59 Epoch: 60 Epoch: 61 Epoch: 62 Epoch: 63 Epoch: 64 Epoch: 65 Epoch: 66 Epoch: 67 Epoch: 68 Epoch: 69 Epoch: 70 Epoch: 71 Epoch: 72 Epoch: 73 Epoch: 74 Epoch: 75 Epoch: 76 Epoch: 77 Epoch: 78 Epoch: 79 Epoch: 80 Epoch: 81 Epoch: 82 Epoch: 83 Epoch: 84 Epoch: 85 Epoch: 86 Epoch: 87 Epoch: 88 Epoch: 89 Epoch: 90 Epoch: 91 Epoch: 92 Epoch: 93 Epoch: 94 Epoch: 95 Epoch: 96 Epoch: 97 Epoch: 98 Epoch: 99 Epoch: 100
Trong khi, Semi-Supervised GMM phân cụm rất chính xác, chỉ nhờ vào 5 điểm cho trước thuộc mỗi phân phối Gaussian